#!/usr/bin/python
#coding=utf-8

from __future__ import division, print_function
from ROOT import TCanvas, TH2D,TH1D, TFile, TNtuple, TF1, TTree, TBrowser, TSpectrum, TGraph
from array import array

##---Calibrazione dei TAC
detectors_names = ["up", "down", "up_and_down"]
delay_up = range(1, 20, 2)
delay_down =  range(0, 21, 2)
delay_up_and_down = range(1, 20, 1)
delays = [delay_up, delay_down, delay_up_and_down]
file_tac = open('./risultati/calibrazione_TAC.dat', 'w')
line = 'time = canale*c +q\n'
file_tac.write(line)


for i in range(3):
    file_name = "dati/calibrazione_tempi_" + detectors_names[i] + ".root"
    file = TFile(file_name)
    tree = file.Get("pjmca")
    isto = TH1D("isto calibrazione", "Calibrazione tempi " +
            detectors_names[i] + ";[channel];[count]", 2200, 0, 2200)
    isto.SetStats(0)
    tree.Project("isto calibrazione", "ch" + str(i))
    s = TSpectrum()
    s.Search(isto, 1, "", 0.75)
    peaks = [s.GetPositionX()[j] for j in range(s.GetNPeaks())]
    peaks.sort()
    ofile = open("dati_graph_cal_" + detectors_names[i] + ".dat", "w")
    for j in range(s.GetNPeaks()):
        line = '{0}\t{1}\n'.format(peaks[j], delays[i][j])
        ofile.write(line)
    ofile.close()
    canv_cal = TCanvas("canv_cal", "canv_cal")
    calibration_graph = TGraph("dati_graph_cal_" + detectors_names[i] + ".dat")
    calibration_graph.SetMarkerStyle(3)
    calibration_graph.SetTitle('Calibrazione TAC {0};[canale];Ritardo [ns]'.format(detectors_names[i]))
    calibration_graph.Draw("ap")
    f = TF1("fit calibrazione", "pol1", 0, 2000)
    canv_cal.SaveAs('./risultati/cal_TAC_{0}.eps'.format(detectors_names[i]))
    fit_par = calibration_graph.Fit(f,'q')
    result = '\nCalibrazione TAC {0}\n'.format(detectors_names[i])
    q = f.GetParameter(0)
    c = f.GetParameter(1)
    chi2 = f.GetChisquare()
    result += 'c = {0}\nq = {1}\nchi2 = {2}\n'.format(c, q, chi2)
    file_tac.write(result)
file_tac.close()
c_up_down = c
q_up_down = q

##---fine Calibrazione dei TAC

##---Calibrazione degli scintillatori
file_name = 'dati/ampia_statistica_con_tracciatori.root'
file = TFile(file_name)
tree_raw = file.Get('pjmca')

isto_up_down_raw = TH2D('isto_up_down_raw', 'Up vs down non calibrato;[canale];[canale]', 2000, 30, 2030, 2000, 30, 2030)
cut_calibrazione = 'ch3>50 && ((ch0<1000 && ch1<850) || (ch0>1200 && ch1>1300))'
tree_raw.Project('isto_up_down_raw', 'ch0:ch1', cut_calibrazione)
canvas_up_down_raw = TCanvas('canvas_up_down_raw', 'canvas_up_down_raw')
isto_up_down_raw.SetStats(0)
isto_up_down_raw.Draw()
canvas_up_down_raw.SaveAs('./risultati/up_vs_down_raw.eps')

#Distribuzione teorica (Tutte le quantità sono in cm)
#Distanza del punto medio dalla faccia dello scintillatore
h = 25
#Distanza del limite dello scintillatore dalla verticale dei tracciatori
d = 10 
#Semi-larghezza dello scintillatore
L = 10.5
#Posizione dei tracciatori
tr1 = 10
tr2 = 110

##---Calibrazione rivelatore UP
#Distanza dallo scintillatore dal punto medio del segmento congiungente le
#facce dei tracciatori vicine tra loro
H = 185
#Raggio del cerchio proiettato dal cono dei tracciatori sullo
#scintillatore
D = 5/h*H

isto_raw_up = TH1D('calibrazione_x_up', 'Calibrazione x, rivelatore UP; [canale]; [count]', 50, 0, 2000)
tree_raw.Project('calibrazione_x_up', 'ch0', cut_calibrazione)

#Definizione della formula per il fit 
x = '([2]*x - [1])'
formula = '[0]/sqrt({0}**2 + {1}**2)*atan({2}/sqrt({0}**2 + {1}**2))'.format(x, H, L)

f1 = TF1('dist_teo',  formula, -d , D)
f1.SetParameters(100000,+400,1)
fit = isto_raw_up.Fit(f1,'RSQ', ' ', 400, 800)

f2 = TF1('dist_teo',  formula, -d , D)
f2.SetParameters(100000,+400,1)
fit = isto_raw_up.Fit(f2,'RSQ+', ' ', 1360, 1760)

x1 = f1.GetMaximumX(200, 700)
x2 = f2.GetMaximumX(1700, 1800, 10)

#La posizione lungo lo scintillatore si trova con x(cm) = c * x(bin) + q
c_up = (tr2-tr1)/(x2-x1)
q_up = -x1*c_up + 10

canvas_raw_up = TCanvas('canvas_raw_up', 'canvas_raw_up')
isto_raw_up.SetStats(0)
isto_raw_up.Draw()
canvas_raw_up.SaveAs('./risultati/cal_x_up.eps')

##---fine Calibrazione rivelatore UP---

##--- Calibrazione rivelatore DOWN
#Distanza dallo scintillatore dal punto medio del segmento congiungente le
#facce dei tracciatori vicine tra loro
H = 85
#Raggio del cerchio proiettato dal cono dei tracciatori sullo
#scintillatore
D = 5/h*H

isto_raw_down = TH1D('calibrazione_x_down', 'calibrazione x, rivelatore DOWN; [canale]; [count]', 50, 0, 2000)
tree_raw.Project('calibrazione_x_down', 'ch1', cut_calibrazione)

#Definizione della formula per il fit 
x = '([2]*x - [1])'
formula = '[0]/sqrt({0}**2 + {1}**2)*atan({2}/sqrt({0}**2 + {1}**2))'.format(x, H, L)

f1 = TF1('dist_teo',  formula, -d , D)
f1.SetParameters(100000, 650, 1)
fit = isto_raw_down.Fit(f1,'RSQ', ' ', 500, 800)

f2 = TF1('dist_teo',  formula, -d , D)
f2.SetParameters(100000, 1700, 1)
fit = isto_raw_down.Fit(f2,'RSQ+', ' ', 1350, 1660)

x1 = f1.GetMaximumX(200, 700)
x2 = f2.GetMaximumX(1500, 1600)

#La posizione lungo lo scintillatore si trova con x(cm) = c * x(bin) + q
c_down = (tr2-tr1)/(x2-x1)
q_down = -x1*c_down + 10

canvas_raw_down = TCanvas('canvas_raw_down', 'canvas_raw_down')
isto_raw_down.SetStats(0)
isto_raw_down.Draw()
canvas_raw_down.SaveAs('./risultati/cal_x_down.eps')

#---fine Calibrazione rivelatore DOWN---

#---fine Calibrazione degli scintillatori---

##---Creazione dell'albero calibrato nel file 'eventi_calibrati.root'
file_calibrato = TFile('eventi_calibrati.root', 'recreate')
tree = TTree('eventi_calibrati', 'Albero con eventi calibrati')
cm_up  = array('d', [0])
cm_down = array('d', [0])       
time_up_down = array('d', [0])       
tracciatori = array('i', [0])       
tree.Branch('x_up', cm_up, 'x_up/D')
tree.Branch('x_down', cm_down, 'x_down/D')
tree.Branch('time', time_up_down, 'time/D')
tree.Branch('trac', tracciatori, 'trac/B')
for i in range(tree_raw.GetEntries()):
    tree_raw.GetEntry(i)
    canale_up = tree_raw.ch0
    canale_down = tree_raw.ch1
    canale_up_down = tree_raw.ch2
    if tree_raw.ch3 < 50:
        tracciatori[0] = 0
    else:
        tracciatori[0] = 1
    cm_up[0] = c_up*canale_up + q_up 
    cm_down[0] = c_down*canale_down + q_down 
    time_up_down[0] = c_up_down*canale_up_down + q_up_down
    tree.Fill()
file_calibrato.Write()

#Istogrammi di posizione calibrati
canvas_up = TCanvas('canvas_up', 'canvas_up') 
isto_up = TH1D('posizione_up', 'Posizione up; [cm]; [count]', 160, -20, 140)
tree.Project('posizione_up', 'x_up')
isto_up.Draw()
canvas_up.SaveAs('./risultati/up_calibrato.eps')
    
canvas_down = TCanvas('canvas_down', 'canvas_down') 
isto_down = TH1D('posizione_down', 'Posizione down; [cm]; [count]', 160, -20, 140)
tree.Project('posizione_down', 'x_down')
isto_down.Draw()
canvas_down.SaveAs('./risultati/down_calibrato.eps')

##Distribuzione della posizione down per gli eventi con posizine up tra 5-10 cm
#isto_pos_down = TH1D('pos_down', 'Posizione down fissata la up; [cm]; [count]',80, -20, 140)
#tree.Project('pos_down', 'x_down', 'x_up<5 && x_up>0')
#canvas = TCanvas('canvas', 'canvas')
#isto_pos_down.Draw()
input()
file_calibrato.Close()
##---fine Creazione dell'albero calibrato nel file 'eventi_calibrati.root'
